Number of relevant directions in Principal Component Analysis and Wishart random matrices. 
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We compute analytically, for large TV, the probability 7(N+ , TV) that a TV x TV Wishart random matrix 
has TV+ eigenvalues exceeding a threshold TV£, including its large deviation tails. This probability plays a 
benchmark role when performing the Principal Component Analysis of a large empirical dataset. We find 
that J'(TV+,TV) « exp(— / 9TV 2 ^(TV+/TV)), where /3 is the Dyson index of the ensemble and t^c( K ) ls a rate 
function that we compute explicitly in the full range < n < 1 and for any f , The rate function ip((n) displays 
a quadratic behavior modulated by a logarithmic singularity close to its minimum K*(Q. This is shown to be a 
consequence of a phase transition in an associated Coulomb gas problem. The variance A (TV) of the number of 
relevant components is also shown to grow universally (independent of £) as A(TV) ~ In TV for large 

TV. 
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Large datasets, such as financial or climate data, typically 
contain a mixture of relevant fluctuations responsible for hid- 
den inherent patterns and irrelevant minor fluctuations. Com- 
pressing the data by filtering out the irrelevant fluctuations, 
while retaining the relevant ones, is crucial for many prac- 
tical applications. The most commonly used technique for 
such data compression is the "Principal Component Analy- 
sis" (PCA) with an impressive list of applications including 
image processingJj]-[3l], biological microarrays J3,|5t], popula- 
tion genetics JauDi finance JalU], meteorology and oceanog- 
raphy [ 10]. The basic idea of PCA is rather simple. Consider 
for instance an experiment where a set of TV observables is 
recorded M times. The collected data can be arranged into a 
rectangular M x N matrix X and adjusted to have zero mean. 
For example, (X)y might represent examination marks of the 
2-th student (1 < i < M) in the j-th subject (physics, math- 
ematics, chemistry, etc.), or the recorded temperatures of the 
i-th largest American city on the j-th measurement day. The 
product symmetric (TV x TV) square matrix C = X T X rep- 
resents the empirical covariance (unnormalized) matrix of the 
data that encodes all correlations. 

In PCA, one collects eigenvalues and eigenvectors of C. 
The eigenvector {principal component) associated with the 
largest eigenvalue Ai gives the direction along which the data 
are maximally scattered and correlated. The scatter and corre- 
lations progressively reduce as one considers lower and lower 
eigenvalues. Next one retains only the top TV + eigenvalues 
(out of TV) and their corresponding eigenvectors, while dis- 
carding all the lower eigenvalues and eigenvectors. Subse- 
quently, the data can be re-expressed in a new basis with the 
following desirable features (i) the dimension of the new ba- 
sis is smaller (i.e. redundant information has been wiped 
out), and (ii) the total variance of data captured along the new 
eigendirections is larger (i.e. only the most important patterns 
present in the initial data have been retained). 

While this procedure is well defined in principle, one imme- 
diately faces a practical problem, namely: what is the optimal 
number of eigenvalues TV + to retain in a given data set? Un- 
fortunately, a sound prescription for the choice of TV + is gen- 



erally lacking, and one has to resort to any of the ad hoc stop- 
ping criteria available in the literature 111 ill . For instance, in 
the commonly used Kaiser-Guttman type rule lfl2ll one keeps 
the eigenvalues greater than a threshold value (N, where the 
threshold may represent, e.g, the average variance of the em- 
pirical dataset. 

However, this simple operational rule has attracted some 
criticism 1 1311 since randomly generated data (with no under- 
lying pattern whatsoever) may also produce eigenvalues ex- 
ceeding any empirical threshold, just by statistical fluctua- 
tions. It is thus compelling to define a statistical criterion able 
to discriminate whether the observed number of eigenvalues 
TV + exceeding a chosen threshold £TV retain inherent patterns 
of the data or just reflect random noise. 

A Bayesian approach provides a possible answer to this 
problem 11411 . Suppose we observe a number TV + of eigenval- 
ues of C exceeding a fixed threshold decided by some stop- 
ping rule. We need to estimate the probability P(Mi\N + ) that 
the data are extracted from a certain model M* (out of a set of 
possible {Mj}), given the observed event TV + . According to 
Bayes' rule: 

P(TV+|M)P(M) 



P(Af,|TV 4 



£,P(JV-+|Af i )P(M i 



(1) 



where P(Mi) is the prior, i.e. the probability that the data 
come from model Mj without the knowledge of the observed 
event TV + , while P(N + \Mi) is called the likelihood of draw- 
ing the value TV + if the data are taken from the model Mj. 
Assuming that one has some apriori knowledge of the priors, 
the most crucial ingredient in Bayesian analysis is the esti- 
mate of the likelihood P(TV+|Mj) for a model Mj. The most 
natural 'null' model to compare data with, in our context, is 
the random model where the 'data' Xij are replaced by pure 
noise, say from a Gaussian distribution with zero mean and 
unit variance. In that case, the corresponding covariance ma- 
trix W = C = X+X is called the Wishart matrix 0]. Thus, 
the natural and important question that one faces is: what is 
the likelihood P(TV + |W) of TV + associated with the Wishart 
matrix model? 
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In this Letter, by suitably adapting a Coulomb gas method 
recently used IU61I17I1 to compute the distribution of the num- 
ber of positive eigenvalues of Gaussian random matrices, we 
are able to compute analytically this likelihood P(N + 1 W) for 
arbitrary £ and large N. To make the N dependence explicit, 
henceforth we use the notation P(N + \W) = T(N + ,N) 
where 7( N + , N) then denotes the full probability distribution 
of the number of eigenvalues of an (JV x N) Wishart matrix 
exceeding a fixed threshold (N. In addition to serving as a 
crucial benchmark for the Bayesian analysis of the number of 
retained components in PCA, we show that the distribution 
y(iV+, N) also has a very rich and beautiful structure. 

Let us first summarize our main results. We consider 
Wishart matrices W = X^X where X is in general a rect- 
angular M x N matrix (M > N) with independent entries 
(real, complex or quaternions, labelled by the Dyson index 
P = 1,2,4) drawn from a standard Gaussian distribution of 
zero mean and unit variance. We focus, for convenience, on 
the case when M — N ~ 0(1)- The matrix W has N non- 
negative random eigenvalues. We compute the distribution 
y(N + , N) for large N where N + is the number of eigenval- 
ues of W exceeding the threshold £N, Setting N + = kN, 
we show that "P(N + = kN, N) behaves for large N as HH 



y(N+ = kN, N) w exp [-/3N 2 ^ c (k)] , < k < 1 (2) 

where the rate function iPq(k) (independent of /3) is com- 
puted analytically for arbitrary £ (see Eq. (fT8l l) and plot- 
ted in Fig. |2] for ( = 1. The rate function has a mini- 
mum (zero) at the critical value n = k*(() where K*(Q — 
(4 arcsec(2/V?) - VC( 4 - C))/2tt. Thus the distribution is 
peaked around N + = k*(£)N which is precisely its mean 
value (N+) — k*(C)N. Around this critical value, the rate 
function is non-analytic and displays a quadratic behavior 
modulated by a universal (^-independent) logarithmic singu- 
larity 



^(0+S)^~(n 2 /2)5 2 /\og\6\ 



(3) 



The physical origin of this non-analytic behavior is traced 
back to a phase transition in the associated Coulomb gas 
problem when k crosses the critical value k*(C). In- 
serting this expression in (ffjl, one finds that, for small 
fluctuations of N + around its mean (N + ) on a scale 
\/ln N, the distribution has a Gaussian form, 7(N + , N) ~ 



exp 



(N+~k*(0N) 2 /2A(N) 




FIG. 1 : Equilibrium density of the Coulomb fluid, Eq. U5\ for C = 1 
(black solid line). The density in general has a disconnected support 
and exists in two different phases, depending on whether k is smaller 
(top panel) or larger (bottom panel) than k* (£). Red dots correspond 
to Montecarlo simulations of a Coulomb gas of 20 charges, where 7 
charges (top) and 12 charges (bottom) were constrained to the right 
°f C = 1 12611 . The corresponding values of the edge point a are 
a = 1.04815 > C (top) and a = 0.701721 < ( (bottom). The area 
under the leftmost blob is 1 — k and under the rightmost blob is n. 



We start by recalling some well known spectral properties 
of Wishart matrices (also known as Laguerre or chiral ensem- 
ble) defined earlier. The N non-negative eigenvalues {A^jof 
the Wishart matrix W are distributed via the joint probability 
density law IU9I1 



?(Ai 



■An) = J-e" 



N 



|A(A)|" (5) 



where A(A) = U j<k ( x j ~ x k\ a = (1 + M — N) — 2//3 
and Zn is the normalization constant. It is also known that 
for large N, Z N ~ exp[~/3n N 2 } where il = 3/4 HEU. 
This joint law (0 can be conveniently recast in the Boltzmann 
form 



where the variance 



T(A 1 ,...,A A r)«ex P (-(/3/2)£;[{A}]) 
where the energy of a configuration {A} is 



(6) 



A(N)=N 2 (( K - K *(C)f) 



fin'- 



■lnJV + 0(l) (4) 
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This leading behavior is also universal, i.e., independent of 
C and is the same as the variance of the number of positive 
eigenvalues in the Gaussian ensembles EH. We thus con- 
clude that on a scale - O(VlniV), the distribution V(N + ,N) 
is Gaussian with mean N + = n*(() and variance in (|4j, but 
with non-Gaussian large deviation tails that are described by 
the general formulae © and cfT~8T >. 



i—l i—1 j^k 

Here, (3/2 stands for the inverse temperature and for simplic- 
ity we focus here on almost-square matrices with M — N 
(and thus a) of 0(1) for large N. This thermodynamical 
analogy, originally due to Dyson 02211 . allows to treat the sys- 
tem of eigenvalues as a Coulomb gas, i.e., a fluid of charged 
particles confined to the positive half-line by two competing 
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interactions: the external linear + logarithmic potential (the 
first two terms) in © tends to push the charges towards the 
origin, while the third term representing mutual logarithmic 
repulsion between any pair of charges tends to spread them 
apart. It is useful first to estimate the typical scale of an 
eigenvalue A for large N. The first two terms typically scale 
for large N as A typ A^ and ln(A typ ) N. In contrast, the pair- 
wise Coulomb repulsion (the third term) typically scales as 
N 2 for large N. Balancing the first and the third term in- 
dicates A typ ~ N for large N. Consequently, as long as 
a ~ 0(1), the second term becomes smaller (~ 0(N)) com- 
pared to the other two terms (~ 0(N 2 )) and hence can be 
neglected. One then expects the average spectral density (nor- 
malized to unity) p(X) = N~ 1 '^ /i (S(X — Aj)) to have the 
scaling form for large N: p(X) — N^ 1 f mp (X/N). The scal- 
ing function / mp (x) = (2n)~ 1 y/ (4 — x)/x is the celebrated 
Marcenko-Pastur (MP) distribution on the compact support 
X 6(0,4)10]. 

Given the joint distribution in (0, our goal is to compute the 
statistics of N+ denoting the number of eigenvalues greater 
than a fixed threshold value QN, i.e., N+ = J2? =1 <?(A 4 ~ (N) 
where 6(x] is the Heaviside step function. The average value 
of N + can be easily estimated from the MP spectral den- 
sity, (N+) = N f*f mp {x)dx = k*{QN where «*(£) = 



(4 arcsec(2/VC) - a/C(4 - C))/2tt. The full probability dis- 
tribution of N+ can be written as the TV-fold integral 

J i=i V »=i / 

(8) 

which we evaluate next for large N by extending the Coulomb 
gas analogy mentioned above. 

The evaluation of the iV-fold integral ([H) in the large 
N limit consists of the following steps: first, we intro- 
duce an integral representation for the S function, 8{x) — 
(2-k)~ x (PN/2) J dpexp(i(3Npx/2). Then, we cast again the 
integrand in the Boltzmann form exp [— (f3/2)E K ({\})} with 



E K ({\}) = E({\}) + AiN (£i S{\ ~ (N) - kN), where 
Ai = \p can be interpreted as a Lagrange multiplier. Written 
in this form, the integral has again a natural thermodynamical 
interpretation as the canonical partition function of a Coulomb 
gas in equilibrium at inverse temperature /3/2. This time, 
however, in addition to the linear confinement and the log- 
arithmic repulsion, the fluid particles (eigenvalues) are sub- 
jected to another (discontinuous) external potential. This ex- 
ternal term involving A\ in the new energy function E K ({\}) 
has the effect of constraining a fraction k of the fluid charges 
to the right of the point x = (N. 

In the large N limit, the Coulomb gas with N discrete 
charges becomes a continuous gas which can be described by 
a continuum (normalized to unity) density function g(X) = 
N- 1 J2lLi S(X-Xi). Consequently, one can replace the orig- 
inal multidimensional integral in ([H} by a functional integral 
over the space of g(X). This procedure, originally introduced 
by Dyson i22ll . has recently been used successfully in a num- 
ber of different contexts (see E [B El H El and refer- 
ences therein). Noting further that the density g(X) is ex- 
pected to have the scaling form g(X) = N^ 1 f K (X/N) for 
large N, each term in the energy E K ({\}) is of the same or- 
der ~ 0(7V 2 ) and is expressed as an integral over the function 
f K (x) (see eq. (TToT > below). The scaled density f K (x) satisfies 
the overall normalization condition J Q dxf K (x) — 1 and the 

constraint dxf K (x) = 1 — k arising from the theta-function 
termini; K ({A}). 

The probability ^ can now be rewritten as y(N + = 
kN,N) = Z K (N)/Z N , where the numerator Z K (N) is the 
following functional integral over f K (x), supplemented by 
two additional integrals over auxiliary variables A\ and A2 
enforcing the two constraints mentioned above: 



Z I; [\)= / ./.I, d.l, (.<•]>]> <J ~N 2 S[f K (x)} 



The action S[f K (x)] is given by: 



(9) 



OO p OO 



§[f K (x)]= / dxxf K (x)- / dxdx' f K (x)f K (x')\n\x-x'\+Ai 
Jo Jo Jo 

I 

Eq. (O is indeed amenable to a saddle point analysis: 



Z K (N) w exp ( -!?-N 2 S[t K (x)} ) , Z N w cxp(~ f3n N 2 ) 



where Qq = 3/4 and f*(x) is obtained from Ss L£afeiJ 



dxf K [x) - (1 - k) \ +A 2 



dxf K {x) - lj 

001 



This gives 

x + A x e{C -x)+A 2 =2 J dx'f*(x') log I a; - x'\ (12) 

for all x € supp(/*(x)) where supp denotes the region when 
charges exist, i.e., f*(x) is positive. Varying the action with 
(11) respect to A\ and A 2 just reproduces the two constraints. The 
= 0. function f*(x) obtained in this way is the rescaled equilib- 
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FIG. 2: The rate function ip^ as a function of ft € [0, 1] from Eq. 
J 1 8t for £ = 1. In the two insets, a sketch of the corresponding 
Coulomb gas equilibrium densities (Eq. J15t ) in the two regimes 
k < (left) and k > (right). Only at k = «*(1) the two 

disconnected blobs of particles on both sides of the barrier at £ = 1 
merge into a single continuous density (the usual Marcenko-Pastur 
distribution / mp (a;)). 



rium density of a Coulomb fluid where a fraction k of parti- 
cles are constrained to lie to the right of a barrier at x = £. It 
is easy to show [26] that f*(x) reduces to the unconstrained 
MP density / mp (a;) when n = «*(£). 

Thus, for large N, the saddle point analysis precisely pre- 
dicts the result in (O with the rate function given by 



(13) 



It then remains to solve the saddle-point equation (fT2l . Taking 
one more derivative with respect to x (x ^ 0): 



i = Pr 
2 



dx' 



(14) 



where Pr denotes Cauchy's principal value, supplemented 
with the constraints J* °° dxf*(x) — 1 and J Q dxf*(x) — 

1 - K. 

To invert the singular integral equation (TT4b is a nontrivial 
challenge. One often needs to guess the solution and verify 
it a posteriori. To get a feeling how the solution may look 
like, we first did a Monte Carlo simulation of the Coulomb 
gas which brought out the following very interesting features. 
For k ^ k*((),0 and 1, the equilibrium charge density f*(x) 

(i) generally consists of two disconnected supports: a blob of 
(1 — k)N eigenvalues to the left of £ separated from a second 
blob of kN eigenvalues to the right of £ (see Fig. [TJ, and 

(ii) the actual shapes of the two blobs depend on whether k < 

(top panel of Fig. [TJ or n > «*(C) (bottom panel of Fig. 
[TJ, i.e. the fluid undergoes a phase transition as k is varied 
across the critical point «*(£)■ When k — > k*(Q, the two 
blobs merge into a single support solution which is precisely 
the MP spectral density f mp (x). 



Extracting this two-support solution for a generic k ^ 
0, 1 from the singular integral equation (Q3J poses the 
main technical challenge that we have succeeded in solving. 
There exists a well known Riemann-Hilbert method B27I1 for 
solving such singular integral equations when the solution 
has a single support 02811 . For two disjoint supports, this 
method ||27|] was recently generalized in a number of different 
problems Ilia Il7l I29W31I1 . Adapting this generalized method 
to our problem (for details see 112611 ). we find that the two- 
support solution of (Q3J satisfying the two constraints is given 
explicitly by the expression 



(a; - a) (4 + C- a-x) 



2tt 



x(x - C) 



(15) 



The edge point a, depending on k, is to be determined as 
the unique solution of the constraint dxf*(x) = 1 — n, 
which can be rewritten in terms of generalized hypergeomet- 
ric functions lf26Tl . We have a > £ for k < K*(C)> a < C 
for k > K*(Q and a = ( for k = K*(£). The density 
f£(x) is supported on the union of two disconnected inter- 
vals for k ^ K*(C), [0, C] U [a, 4 + ( - a] (if k < K*(Q) or 
[0, a] U [C, 4 + C - a] (if k > n*(Q) [see, e.g., the insets of 
Fig. [2). For k = k*(C)^ the two intervals merge into a single 
one and j*{x) = f mp (x) as expected. In the extreme cases 
k = 0,1, the density has also a single support and its ana- 
lytical form is known from earlier works 112 111 . The function 
/*(x) is plotted in Fig. [ljfor a specific choice £ = 1, together 
with the result of Montecarlo simulations of the Coulomb gas, 
showing excellent agreement. 

Next we substitute this explicit solution f*(x) in (QJjJ to 
evaluate the saddle-point action §[f*(x)}. For this, we need 
to compute the single and double integrals in ( TTOj . The dou- 
ble integral can be written in terms of a simple integral after 
multiplying (Q~2j by f*(x) and integrating over x. This leads 
to: 



Sffl = ^i(0-4(1-K)-A a ] 



(16) 



where /-t n (C) = Jo°° a:; ™/£( a; )^a ; is the n-th moment of the 
density. The first moment can be explicitly computed as 

/Xi(C) = 1 - a + x + C _ S l- It remains to fix the two 
Lagrange multipliers A\ and Ao, assigning suitable x values 
in Eq. (Q~2j. Skipping details B26I1 . it turns out that A\ and A 2 
can be expressed in terms of a pair of functions: 



- ± - 



1 (a-x)(4 + (-a-x) 



x(x - C) 



(17) 



defined for x ^ supp(/*(a;)). It can be shown that these func- 
tions are essentially the moment generating functions of /* (x) 

Combining dT6l > and ( TOJ . the final expression for the rate 
function ip((K) reads: 
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HO - (1 - «) (a - C - 2 £ F (+) (z)^ 
A(C) - (1 - «) (C- a -2^^) 



2 ln(4 + C - a) + 2 j£ c _ a ^ (x) - ±) 
21n(4 + C-a) + 2/~ c _ a d a! (F ( _ ) ( a: )-i) 

I 



K < K*(C) 

08) 



where /2(C) = 2 - 2a + a 2 /4 + 2( - </4. The behavior of 
tp£ (k) around the minimum n = n* (£) in (01 is found after a 
lengthy but straightforward expansion i26ll . 

In summary, using a Coulomb gas approach we computed 
analytically the likelihood P(iV + |W) of retaining N + princi- 
pal components if a completely random (pure noise) N x N 
prior W (Wishart random matrix) is assumed for the underly- 
ing data. This likelihood is an essential ingredient in Bayes' 
formula (|T), which in turn gives the probability that a certain 
number N + of significant eigenvalues is observed if the data 
represents pure noise. Thus the rate function computed in ( [T8T l 
serves indeed as a benchmark to gauge the significance of re- 
sults obtained for empirical data. We also found an interesting 
phase transition in the associated Coloumb gas problem when 
two separated blobs of Coulomb charges merge onto a single 
one as a critical value N + /N = n = k*(C) is approached. 
This phase transition is responsible for the appearance of a 
logarithmic singularity in the rate function close to its min- 
imum, which in turn leads to a highly universal logarithmic 
growth of the variance of N+ with N. 

The present work can be developed in several directions: on 
one hand, it would be interesting to analyze the case N ^ M 
and check if the universality found for N ~ M still persists. 
On the other hand, improved stopping criteria based on null 
random matrix results are very much called for. The Coloumb 
gas method developed here for null random matrices may in- 
deed be useful in that direction. 

S.N.M. acknowledges the support of ANR grant 2011- 
BS04-013-01 WALKMAT. 



[1] S. S. Wilks, Mathematical Statistics (John Wiley & Sons, New 
York, 1962). 

[2] K. Fukunaga, Introduction to Statistical Pattern Recognition 

(Elsevier, New York, 1990). 
[3] L. I. Smith, A tutorial on Principal Components Analysis 

(2002). 

[4] N. Holter et al., Proc. Natl. Acad. Sci. USA 97, 8409 (2000). 
[5] O. Alter et al., Proc. Natl. Acad. Sci. USA 97, 10101 (2000). 
[6] L.-L. Cavalli-Sforza, P. Menozzi, and A. Piazza, The History 



and Geography of Human Genes (Princeton Univ. Press, 1994). 

[7] J. Novembre and M. Stephens, Nature Genetics 40, 646 (2008). 

[8] J.-P. Bouchaud and M. Potters, Theory of Financial Risks (Cam- 
bridge University Press, Cambridge, 2001). 

[9] Z. Burda and J. Jurkiewicz, Physica A 344, 67 (2004); Z. Burda, 
J. Jurkiewicz, and B. Waclaw, Acta Phys. Pol. B 36, 2641 

(2005) and references therein. 

[10] R. W. Preisendorfer, Principal Component Analysis in Meteo- 
rology and Oceanography (Elsevier, New York, 1988). 

[11] D. A. Jackson, Ecology 74, 2204 (1993). 

[12] L. Guttman, Psychometrika 19, 149 (1954); N. Cliff, Psycho- 
logical Bulletin 103, 276 (1988). 

[13] G. D. Grossman and D. M. Nickerson, Ecology 72, 341 (1991). 

[14] J. O. Berger, Statistical Decision Theory and Bayesian Analysis 
(Springer Series in Statistics, 2nd ed., Springer- Verlag, 1985). 

[15] J. Wishart, Biometrika 20, 32 (1928). 

[16] S. N. Majumdar, C. Nadal, A. Scardicchio, and P. Vivo, Phys. 

Rev. Lett. 103, 220603 (2009). 
[17] S. N. Majumdar, C. Nadal, A. Scardicchio, and P. Vivo, Phys. 

Rev. E 83, 041105 (2011). 
[18] The symbol ~ stands for the precise logarithmic law 

lim Jv ^oolnJ > ( K Ar,Ar)/(/3Af 2 ) = -ip C ( K )- 
[19] A.T. James, Ann. Math. Statistics 35, 475 (1964). 
[20] M. L. Mehta, Random Matrices (Academic Press, Boston, 

1991). 

[21] P. Vivo, S. N. Majumdar, and O. Bohigas, J. Phys. A: Math. 

Theor. 40, 4317 (2007). 
[22] F. J. Dyson, J. Math. Phys. 3, 140 (1962). 
[23] V. A. Marcenko and L. A. Pastur, Math. USSR-Sb 1, 457 

(1967). 

[24] D. S. Dean and S. N. Majumdar, Phys. Rev. Lett. 97, 160201 

(2006) . 

[25] D. S. Dean and S. N. Majumdar, Phys. Rev. E 77, 41 108 (2008). 
[26] See Supplementary Material. 

[27] E. Brezin, C. Itzykson, G. Parisi, and J. B. Zuber, Commun. 
Math. Phys. 59, 35 (1978). 

[28] F. G. Tricomi, Integral Equations (Pure Appl. Math V, Inter- 
science, London, 1957). 

[29] P. Vivo, S. N. Majumdar, and O. Bohigas, Phys. Rev. Lett. 101, 
216809 (2008). 

[30] P. Vivo, S. N. Majumdar, and O. Bohigas, Phys. Rev. B 81, 
104202(2010). 

[31] K. Damle, S. N. Majumdar, V Tripathi, and P. Vivo, Phys. Rev. 
Lett. 107, 177206 (2011). 



